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1. Introduction 

Plain weave composites, reinforced by mutually interlaced systems of unidirectional 
fiber tows bonded to a matrix, belong to a progressive material systems with 
widespread applications in virtually all areas of engineering. As a particular example, 
consider carbon-carbon (C/C) composites, originally developed for the space and 
automobile industry, which now find uses in the medicine owing to their appealing 
biological compatibility with a living soft tissue, e.g. (Pesakova et al. 2003). A proper 
characterization of these material systems, especially from the mechanical response point 
of view, thus appears rather important. 

While a detailed two-dimensional (2D) analysis of a heat conduction problem for 
the evaluation of effective (macroscopic) thermal conductivities seems to be sufficient, 
e.g. (Tomkova 2006, Tomkova et al. 2008), a reliable estimate of the mechanical 
response of such systems requires in general a solution of a full three-dimensional 
(3D) problem. This task, however, presents a significant challenge even if limiting our 
attention to a linear elastic behavior. Not only the characteristic 3D structure of textile 
composites, but also various types of imperfections in woven path developed during the 
manufacturing process preclude a direct formulation of a simple computational model. 

A considerable research effort has been invested in the last two decades into 
providing a simple yet accurate scheme for the predictions of macroscopic elastic 
properties of woven composites. With an increasing level of sophistication, these models 
include modified rule of mixtures, approaches based on classical laminate theories (CLT) 
and detailed three-dimensional finite element method (FEM) based simulations, see 
e.g. (Cox & Flanagan 1997, Chung & Tamma 1999, Takano et al. 1999, Lomov 
et al. 2007) for a review and comparison of individual approaches. The latter class of 
computational models is considered to be the most accurate one particularly if combined 
with concise geometrical data (Barbero et al. 2006, Zeman & Sejnoha 2004, Lomov 
et al. 2007). 

The FEM simulations show, however, certain disadvantages. Perhaps the most 
critical one is a relatively high computational cost due to laborious preparation of 
finite element meshes. Moreover, to incorporate at least the dominant microstructural 
imperfections observed in real systems into the FEM model is far from being 
trivial and deserves a special treatment typically based on an appropriate statistical 
characterization (Zeman & Sejnoha 2004, Zeman & Sejnoha 2007). The CLT approaches 
are, on the other hand, easy to implement and provide a reasonable approximation of the 
in-plane elastic moduli. However, since this class of models approximates the composite 
as a coupling of serial and parallel laminates stacked to resemble the actual geometry, it 
becomes inadequate when predicting the out-of-plane response. Therefore, a procedure 
offering a reasonable compromise between the accuracy of FEM-based modeling and 
simplicity of traditional CLT methods is still on demand. 

In the last decade, effective media theories, widely used in continuum 
micromechanics (Bohm 2005), have been recognized as an attractive alternative to 
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CLT-based methods. Such an approach was pioneered by Gommers et al. (1998) 
and Huysmans et al. (1998), who modeled knitted composites as an assembly of 
spherical fibers in an isotropic matrix and used the Mori-Tanaka (M-T) method (Mori 
& Tanaka 1973) to evaluate the overall response. Further advancements in the 
field include the Transformation Field Analysis approach due to Bahei-El-Din et al. 
(2004) and the work of Barbero et al. (2005) in the framework of the theory of 
periodic eigenstrains. All these studies report good correspondence with experimental 
data with an error comparable to experimental scatter. Moreover, with regard to 
imperfect textiles, the Mori-Tanaka method appears particularly useful as it allows, 
through the application of orientation averaging techniques, see e.g (Yushanov & 
Bogdanovich 1998, Gommers et al. 1998, Schj0dt Thomsen & Pyrz 2001, Duschlbauer 
et al. 2003, Jing et al. 2003, Doghri & Tinel 2006, and references therein), for a direct 
introduction of imperfections in the fiber-tow path represented here by histograms of 
distribution of the fiber-tow orientation angles. It is worth noting that such histograms, 
when constructed for all plies in the laminate, also reflect, at least to some extent, tow 
path imperfections due to inter-layer shift typical for real material systems displayed in 
Figure l(a),(b). 

A successful application of the M-T method for the prediction of the effective 
material parameters of textile composites including the above knowledge of the actual 
microstructure requires, however, completion of the following tasks: 

• Quantification of the real micro (meso) structure through a detailed evaluation of 
images of real material samples. This part of the analysis is briefly addressed in 
Section 2 for a C/C composite specimen. 

• The basic geometrical information are then used to construct an idealized three- 
dimensional periodic unit cell exploiting the geometrical model proposed by Kuhn 
& Charalambides (1999). Such a unit cell serves as a point of departure for a 
subsequent application of the M-T method. Review of the model together with 
essential steps of the FEM based simulations using the first-order homogenization 
technique is provided in Section 3. 

• Formulation of the Mori-Tanaka method is then presented in Section 4 with 
emphasis given to the symmetry of the overall material stiffness matrix and 
capturing interaction between individual tows. It is shown that special care is 
required when replacing the actual fiber tow by an equivalent ellipsoidal inclusion. 
In the present formulation, the shape of the equivalent ellipsoid is thus treated as 
an internal parameter of the method determined by matching the M-T estimates 
with the results derived in Section 3. 

• Two possible approaches to the calibration of the internal parameters of the method 
are considered. In the first variant, the optimal ellipsoidal shape is found by 
matching directly the results of the finite element simulations, executed on a 
"training" set reflecting the in-situ determined scatter of geometrical parameters. 
It is worth noting that such procedure allows us to introduce possible deviations 
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from the "ideal" (average) geometry of the textile structure presented in Section 2. 
The second approach employs the finite element data at hand to propose a simple 
relation between basic parameters of the textile composite and parameters of the 
optimal ellipsoid. Section 5 is concluded by verification and validation of the 
developed heuristics. 

• Once the shape of the ellipsoid is calibrated, the orientation averaging can be used 
in conjunction with the histograms of the fiber-tow orientation angle. This final 
step leading to estimates of the overall elastic response of imperfect C/C textile 
composites is examined in Section 6. Section 7 then summarizes the final results 
and offers possible extensions particularly with account to intrinsic porosity of these 
material systems. 

In the following text, the Voigt representation of symmetric tensorial quantities 
is systematically employed, e.g. (Bittnar & Sejnoha 1996). In particular, a, a and 
A denote a scalar value, a vector or a matrix representation of a second-order tensor 
and a matrix representation of a fourth-order tensor, respectively. Other symbols and 
abbreviations are introduced in the text as needed. 

2. Microstructure evaluation 

As already mentioned in the introductory part, obtaining reliable predictions of the 
effective mechanical properties of textile composites requires a thorough analysis of 
their actual microstructure. Figure 1(a) shows a particular C/C composite laminate 
consisting of eight layers of carbon fabric Hexcel G 1169 bonded to a carbon matrix. A 
total of twenty such specimens having dimensions of 25 x 2.5 x 2.5 mm were fixed into 
the epoxy resin and after curing subjected to final surface grounding and polishing using 
standard metallographic techniques to produce specimens suitable for the subsequent 
image analysis. 

While image analysis software LUCIA G® might be used directly to process the 
actual image of the specimen in Figure 1(a), it proves more advantageous, owing to a low 
color contrast of the carbon reinforcement and carbon matrix, to collect the necessary 
geometrical information from its binary counterpart plotted in Figure 1(b). Several such 
sections taken from various locations of the laminated plates were examined to obtain 
basic statistics of various parameters including segment dimensions, fiber tow thickness, 
shape of the fiber tow cross-section, etc. The resulting values are stored in the second 
column of Table 1. The averages of basic geometrical data were finally used to construct 
an equivalent or rather an ideal periodic unit cell (EPUC) appearing in Figure l(c,d) 
employing the description due to Kuhn & Charalambides (1999). The three-dimensional 
geometric model is defined by four parameters: the tow wavelength 2a, the tow height 
6, tow spacing g and the layer thickness h, cf. Figure 1(c). 

Note, however, that the real composite shows a number of imperfections which 
certainly should not be completely disregarded. It will be seen later in Section 4 that the 




Figure 1. Equivalent periodic unit cells; (a) color image of real composite sample, (b) 
binary image, (c) cross-section of an equivalent periodic unit cell, (d) three-dimensional 
view, (e) approximation of centerlines (Vopicka 2004), (f) distribution of inclination 
angles. 



nonuniform waviness and to some extent also the fiber inclinations due to production- 
related mutual shift of individual layers clearly visible in Figure 1(b) can be accounted 
for by utilizing histograms of inclination angles derived from centerlines of individual 



Table 1. Quantification of microstructural parameters 




Carbon/ Carbon 


E-glass/Vinylester 


E-glass/Epoxy 


Parameter 


Tomkova (2004) 


Scida et al. (1999) Kolleg 


;al & Sridharan (2000) 


a [/iin] 


2,250 ± 155 


1,200 


620 


h [/im] 


300 ± 50 


X 


X 


b [/im] 


150 ± 20 


50 


100 


g [fim] 


400 ± 105 


20 


20 


Ctow [%) 


53.2 ± 1.8 


79.8 


69.7 
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fiber tows, see Figure l(e,f) and (Vopicka 2004) for more details. The idealized geometry 
in Figure 1(c) assumes, nevertheless, the centerlines of the warp and fill systems of tows 
in a simple trigonometric form (Kuhn & Charalambides 1999) 



3. Periodic unit cell analysis 

Having quantified the real microstructure, the resulting EPUC can be readily employed 
to provide FEM estimates of the required effective moduli. This particular step of the 
proposed analysis scheme will now be briefly reviewed. 

To that end, consider an EPUC in Figure 1(d) with the local coordinate system 
defined such that the local x\ axis is aligned with the fiber tow direction. Definitely 
the most tedious step in the entire analysis is preparation of a three-dimensional finite 
element mesh complying with the periodic boundary conditions (the same positions 
of the element nodes on the opposite faces of the cell). Here, the elements of CAD 
operations combined with volumetric modeling capabilities of ANSYS® package are 
used to generate the finite element mesh employing the mapped meshing technique 
discussed by Wentorf et al. (1999) and Matous et al. (2007). 

In order to ensure symmetry of the resulting FEM mesh, a primitive block of 
the tow shown in Figure 2(a) is modeled first. Next, using mirroring, copying and 
merging operations, the whole volume of one reinforcement layer is generated. Finally, 
the volume corresponding to the matrix phase is generated by subtracting the body 
of reinforcements from the matrix as depicted in Figure 2(b). To reflect the required 
periodicity only one of the two opposite faces is meshed using the advancing front 




(1) 



(a) 




(c) 




Figure 2. Finite element mesh generation; (a) CAD model of primitive volume, 
(b) CAD model of PUC, (c) FEM mesh of fiber tows, (d) FEM mesh of PUC. 
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technique and then copied to the associated one. At last, the tetrahedral elements 
corresponding to tows and matrix are generated based on the data created in the 
previous steps, leading to finite element meshes shown in Figure 2(c,d). 

The further numerical treatment now proceeds as follows, cf. (Michel et al. 1999). 
Suppose that the periodic unit cell in Figure 2(d) is loaded by a macroscopic strain 
vector E. In view of the assumed microstructure periodicity, the local displacement 
field u then admits the following decomposition 

u(x) = X(x)E + u*(x), (2) 

where u* represents a periodic fluctuation of u due to the presence of heterogeneities 
and matrix X stores the coordinates of x. The local strain then assumes the form 

e(x) = E + e* (x), (3) 

where the fluctuating part e* vanishes upon the volume averaging. Next, introducing 
Equation (3) into the principle of virtual work (the Hill-Mandel lemma) yields 

(5s t (x)(t(x)) = (5e eT (x)(T e (x)) = (5e^ \x)cr\x)) = 0, (4) 

where ( ) stands for the volumetric averaging with respect to the PUC and - £ is used to 
denote a quantity in the local coordinate system. The local stress field then reads 

<T t (x) = L t (x)(E t + e* e (x)), (5) 

where L e is the material stiffness matrix. Relating the strains in the local and global 
coordinate systems by the well-known relations E = T £ E, e l = T £ e, see e.g. (Bittnar 
& Sejnoha 1996), and inserting Equation (5) into Equation (4) yields the stationarity 
conditions in the form 

(5e* T (x)T T £ (x) \L\x)T £ (x) (E + e*(x))]) = 0, (6) 

to be satisfied for all kinematically admissible variations 5e*. 

The homogenized stiffness matrix £ FEM follows from post-processing of the solution 
of six independent elasticity problems, discretized using conforming FEM procedure, 
see (Zeman 2003, Zeman & Sejnoha 2004) for further details. In particular, each column 
of L FEM coincides with the volume averages of local stress cr resulting from a macroscopic 
strain with one component set to one and with the remaining entries equal to zero. 



4. Application of the Mori-Tanaka to woven composites 

4-1. Overall stiffness of composite with non-aligned inclusions 

Consider an A-phase composite with an isotropic matrix phase having the stiffness 
matrix L and being reinforced with (N — 1) families of ellipsoidal heterogeneities. 
Each heterogeneity is characterized by the stiffness matrix L r and occupies a volume 
Q r . With reference to (Benveniste et al. 1991), the Mori-Tanaka estimate of the overall 
stiffness matrix L M ~ T then reads 

L M ~ T = L + (j^ c r (L r - L ) T)j (col + £ c r T^ , (7) 
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where c r denotes the volume fraction of the r-th phase. The corresponding partial strain 
concentration factor T r has the form 

T r = (I + P r (L r -L ))- 1 , (8) 

where the P r matrix is provided by 

P r = J q To(x-x')dx'. (9) 

Function To is related to Green's function of an infinite medium with stiffness matrix 
L (see, e.g. (Ponte Castaneda & Willis 1995, Section 3.1) for more details). It follows 
from the celebrated work of Eshelby (1957) that for ellipsoidal inclusions, P r is constant 
and can be evaluated as 

Pr = S r L~\ (10) 

where S r is the Eshelby matrix. When the matrix phase is isotropic, explicit expressions 
for S r can be found in, e.g. (Eshelby 1957, Mura 1987). 

While the M-T model has proved itself to be accurate for composites reinforced 
either by randomly oriented or aligned inclusions with an identical shape, in general 
case it may lead to a non-symmetric stiffness matrix L M ~ T , see e.g. (Benveniste 
et al. 1991, Ferrari 1991, Ponte Castaneda & Willis 1995) for the in-depth discussion. 
In this work, a simple re-formulation proposed by Schj0dt Thomsen & Pyrz (2001) is 
employed to preserve the overall symmetry of the stiffness matrix using the orientation 
averaging. 

To this end, we approximate the material system under investigation as a two- 
phase composite (N = 2) consisting of an isotropic matrix (r = 0) and with index r = 1 
collectively denoting the reinforcing tow phase, composed of heterogeneities of identical 
shape but different orientations. | Suppose for a moment aligned heterogeneities . Then, 
the overall stiffness matrix is symmetric (Benveniste et al. 1991) and can be decomposed 
to 

L M ~ T = L + Cl ((L, - L ) Ti) ((1 - ci)I + cxTx)- 1 = + Cl L«. (11) 

Notice that due to assumed isotropy of the matrix phase, the matrix is independent 
of the reference coordinate system while stores the orientation-dependent part. 
Following (Schj0dt Thomsen & Pyrz 2001), the estimate of the overall stiffness of a 
composite reinforced with non-aligned heterogeneities is provided by 

^-^I^+dp)), (12) 

where the double brackets (( )) denote averaging over all possible orientations. In 
particular, when the orientation of each heterogeneity is parametrized in terms of the 

| Therefore, the employed geometrical data reduce to the volume fraction of one of the phases and 
appropriate quantification of orientation distribution. 
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Euler angles 0, 9 and (3, see Figure 3,§ the orientation-dependent part can be expressed 
as 

L^(9, 0, (3) = T](6, 0, (3)L^(0, 0, 0)T £ (9, 0, (3), (13) 

where the explicit expression of the transformation matrix can be found in 
e.g. (Schj0dt Thomsen & Pyrz 2001) and (Zeman 2003, Appendix A). 




Figure 3. Definition of the Euler angles. 



The orientation average then follows from 

((L^ )) = fj jT jT (9, 0, (3)g(9, 0, (3) d6 d0 d(3, (14) 

with g(9,cf),(3) denoting the joint probability density describing the distribution of 
individual angles. 



4-2. Application to plain weave composites with ideal geometry 



To examine the theoretical formulation presented in the previous Section, consider 
again an ideal plain weave textile composite already studied in Section 3. In this 
particular case, the joint probability density function g(6,(j),{3) results from the 
harmonic shape of the centerline, recall Equation (1). Applying the change of variable 
formula (Rektorys 1994, Section 33.9), we obtain after some algebra the expression of 
the probability density in the form 



g(9,4>,(3) = { 



2a l + tan 2 (#) 



* ^/bH 2 -4a 2 tan 2 (#) 

otherwise, 



if = 0, (3 = and - a < 9 < a, 



where 



a = arctan 



2a 



§ Note that so-called "x 2 convention" is used; i.e. a conversion into a new coordinates system follows 
three consecutive steps. First, the rotation of angle <j> around the original X% axis is done. Then, the 
rotation of angle 6 around the new x 2 axis is followed by the rotation of angle (3 around the new x 3 
axis to finish the conversion. 
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Equation (14) then becomes 



£(i,war P )\\ = 2«T l + tan 2 (fl) =L m {odo)d ^ (15) 



J -<* yW 2 -4a 2 tan 2 (0) 
and similarly for the fill system we get 

XX 77 77 •'T/a-a ^&% 2 -4a 2 tan 2 (0) 2 

Following Equation (12), the resulting homogenized stiffness matrix of a plain weave 
composite then reads 

L M - T ( Cl ,aA£o,£i,Si) = £ {0) + | (((£ (1 ' fiU) » + ((£ (1 < warp) ))) , (17) 

where the averages of the basic geometrical parameters a, b and the material parameters 
of the two-phase composite including the volume fraction of individual phases are 
assumed to be known quantities. The matrices ^2,( 1 ' warp )^> and (^L^ 1 ' 611 ^ implicitly 
depend, however, on the Eshelby matrix Si, which is yet to be determined. 

To take advantage of the closed-form Eshelby solution (Eshelby 1957), it will be 
assumed that the actual shape of the fiber tow can be well represented by an equivalent 
ellipsoid with semi axes £i > £ 2 > £3 > 0. Then, the accuracy of the M-T method is 
governed by a proper choice of the semi-axes as exemplified by the following case study. 
In particular, three representations is considered: (i) a spherical shape (£1 = £ 2 = £3), 
(ii) a cylindrical shape (£1 — + 00, £2 = £3) and (iii) an ellipsoid (£1 = 1, £2 = 0.5, £3 = 0.1). 
The unit cell with average geometrical parameters appearing in the second column of 
Table 1 and the constituent properties stored in the second column of Table 2, i.e. 
C/C composite system, are considered. Note that in order to achieve the maximum 
phase stiffness contrast, the tow parameters shown in Table 2 correspond to the pure 
carbon fibers. The corresponding homogenized stiffness matrix entries are stored in 
Table 3 together with the FEM data. 



Table 2. Material parameters of individual phases 


in local coordinate system 




Carbon/Carbon 


E-glass / Vinylester 


E-glass/Epoxy 






Barbero et al. (2005) 


Barbero et al. (2005) 


Matrix 








E [GPa] 


30 


3.4 


3.12 


V 


0.19 


0.35 


0.38 


Fiber tow 








E A [GPa] 


210 


58.397 


51.352 


G A [GPa] 


86 


8.465 


5.342 


E T [GPa] 


72 


20.865 


15.040 


G T [GPa] 


27.7 


7.527 


5.342 


V A 


0.27 


0.241 


0.262 
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Table 3. Homogenized stiffness matrix entries (C/C system) 



Method 


L n 


L\2 


Liz 


-^33 


L44 


Lqq 




[GPa] 


[GPa] 


[GPa] 


[GPa] 


[GPa] 


[GPa] 


FEM (Ideal geometry) 


89.94 


16.55 


14.06 


48.67 


20.26 


41.53 


M-T (Spherical inclusion) 


64.73 


14.57 


14.85 


54.56 


24.17 


28.71 


M-T (Cylindrical inclusion) 


103.6 


16.15 


15.55 


52.74 


23.99 


28.75 


M-T (Ellipsoidal inclusion) 


88.35 


16.85 


14.94 


50.24 


21.54 


41.19 



Clearly, comparing the M-T estimates with FEM based results allows us to draw the 
following two conclusions: (i) the Mori-Tanaka method appears as a reliable alternative 
to the first-order periodic homogenization based on the finite element method, (ii) the 
choice of the Eshelby matrix can hardly be made arbitrarily. It should be noted that 
the illustrative results correspond to volume fraction c\ ~ 50%, for which the accuracy 
of the Mori-Tanaka method typically deteriorates. Therefore, the optimal ellipsoid not 
only accounts for the tow geometry, but also for approximately captures tow interactions 
due to non-dilute volume fractions of the reinforcing phase. 

5. Optimal shape of equivalent ellipsoid 

5.1. FEM-based calibration 

The essential goal now becomes to find the optimal shape of the ellipsoid by matching 
the FEM results with the M-T predictions for C/C material system. To take into 
account the observed uncertainties in the textile geometry, a collection of PUCs, rather 
than a single one, is used for the calibration. To generate such a set we exploit the scale- 
invariance of the first order homogenization and set a = 1. The remaining parameters 
were generated using the Latin Hypercube Sampling method (Iman & Conover 1980), 
assuming uniformly distributed random variables with the statistics stored in the second 
column of Table 1. Twenty such unit cells were generated and subject to the FEM-based 
homogenization procedure, yielding the homogenized stiffnesses listed in Table 4. 



Table 4. Summary of homogenized effective properties of training set (C/C system) 



Statistics 


r FEM 


r FEM 
^12 


7- FEM 
^13 


7- FEM 
-^33 


7- FEM 


7- FEM 
-^66 




[GPa] 


[GPa] 


[GPa] 


[GPa] 


[GPa] 


[GPa] 


Average 


86.96 


16.00 


13.65 


47.73 


19.78 


39.13 


Standard deviation 


2.29 


0.40 


0.30 


0.68 


0.34 


1.87 


Optimized M-T (ideal geometry) 


88.81 


16.13 


13.89 


47.35 


20.17 


40.40 


Optimized M-T (training set) 


88.23 


16.78 


14.94 


49.09 


20.09 


40.26 
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Figure 4. Objective function for a single PUC. 



Since the Eshelby matrix depends on the mutual ratio of the ellipsoid semi-axes 
only, it is possible to set = 1, which leaves us with only two parameters undetermined. 
To characterize the discrepancy between the FEM and M-T solution, the following error 
measure is introduced 



max 

i,j=l,...,6 



-FEM 



-L 



M-T 

ij 



E(L FEM ,L M ~ T ) 

When n PUCs are considered, the objective function assumes the form: 



\ 



FEM 
(i) > 



■M-T 



;is) 



(19) 



where the superscript z represents the i-th member in the training set. The optimal 
shape characterized by £| and £3 can be then found from the minimization procedure 

G arg min F(£ a ,6). (20) 

A graphical representation of the objective function assuming a single (average) periodic 
unit cell is plotted in Figure 4 for the sake of illustration. 

Note that the explicit expression of the Eshelby matrix is not available in 
this particular case, which essentially precludes the use of classical gradient-based 
optimization algorithms. The stochastic optimization methods, on the other hand, 
appear to be a more appropriate choice. The particular algorithm, based on the 
surrogate function model combined with evolutionary algorithm adopted in the present 
study is briefly described in Appendix A. 

The solution of the optimization problem then yields the optimal values of semi- 
axes £2 = 0.486 and £3 = 0.092 with the optimal value F* = 1.6. It is worth noting that 
the optimum compares rather well with the case when the minimization is performed 
with respect to the ideal unit cell only, for which E* = 1.3, see also Table 4 for a 
comparison in terms of stiffness matrix entries. This confirms the predictive capabilities 
of the M-T approach, at least in the range of addressed geometrical variations. 
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5.2. Heuristic calibration 

Although the advocated M-T approach seems to offer an efficient way to the prediction 
of the homogenized properties, especially when handling composites with random tow 
imperfections, cf. Section 6, it still requires the reference FEM simulations to tune the 
Eshelby matrix. Therefore, it appears advantageous to establish a heuristic link between 
the EPUC parameters and optimal ellipsoidal shape. In previous works (Gommers 
et al. 1998, Huysmans et al. 1998), such a relation was derived from the local centerline 
curvature and calibrated using selected elastic constants. The current framework, on 
the other hand, offers a possibility to systematically use information contained in the 
previously generated set of EPUCs. 

In the first step of the analysis, the optimization procedure is executed 
independently for each EPUC (i.e. with objective function (18)), yielding a set of 
optimal parameters {£2(1)' £3(1)}^ = 1,2,..., 21. Subjecting the results to correlation 
analysis, see e.g. (Rektorys 1994, Section 34.5), reveals that the parameter is strongly 
correlated with g/a ratio (with the coefficient of correlation equal to —0.8), while it is 
almost independent of b/a value. An analogous trend can be observed between Q and 
b/a parameter. Such results authorize us to postulate a simple linear relation between 
the optimal ellipsoid shape and EPUC parameters. The optimal fit, now determined 
using objective function (20), finally leads to a semi-empirical formula 

,£ 3 *~---. (21) 

No doubt, such heuristics still builds on a representative finite element simulations 
and as such requires to be verified and validated against independent data. In the 
current work, we examine two plain weave composite systems thoroughly analyzed 
in (Barbero et al. 2005, Barbero et al. 2006): (i) E-glass/ Vinylester composite (Scida 
et al. 1999), (ii) E-glass/Epoxy material system (Kollegal & Sridharan 2000). The 
corresponding geometrical data are stored in Table 1, while the material parameters of 
individual constituents are available in Table 2. It is worth noting that the considered 
material systems offer a considerably different tow volume fractions and elastic constants 
of individual phases than in the calibration step. Moreover, to keep the validation 
objective, the comparison will now be based on orthotropic engineering moduli (see, 
e.g. (Bittnar & Sejnoha 1996)) rather than stiffness matrix entries. 

For the E-glass/ Vinylester composite, performance of the M-T method is compared 
with the Periodic Microstructure Model (PMM) (an alternative micromechanics-based 
method based on a detailed geometrical model due to Barbero et al. (2005)), independent 
finite element study in ANSYS® and experimental data. Results of the comparison, 
reported in Table 5, demonstrate a reasonable match between the M-T predictions and 
remaining values. Although the accuracy of the Young moduli is somewhat inferior with 
respect to detailed models, the shear behavior is predicted very well and also the values 
of the Poisson ratios are consistent with the results of alternative numerical approaches. 

Similar conclusions can be made for the E-glass/Epoxy textile system, see Table 5, 
where even closer match between the detailed numerical model can be observed. In 
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overall, the presented data provide an evidence that the heuristic relation (21) leads to 
reasonably accurate estimates of the homogenized elastic properties. 

6. Application to real geometry of C/C composite system 
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Figure 5. Real measured histograms of inclination angles (Vopicka 2004). 



Table 5. Verification and validation of Mori-Tanaka against Periodic media method 
due to Barbero et al. (2005), Finite clement simulation simulation from (Barbero 
ct al. 2006) and experimental data from (Scida et al. 1999, Kollcgal & Sridharan 2000). 

E-glass/Vinylester E-glass/Epoxy 









Experiment 


PMM 


FEM 


M-T 


Experiment 


PMM 


M-T 
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— E 2 2 


[GPa] 


24.8 ± 1.1 


25.1 
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25.8 


19.29 


18.9 
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E33 


[GPa] 




8.5 ± 2.6 
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12.4 
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8.74 


8.83 


G23 


= G13 


[GPa] 


4.2 ± 0.7 


2.91 


3.16 


4.08 


X 


2.57 


2.92 


G12 


[GPa] 




6.5 ± 0.8 


4.37 


5.52 


6.44 


3.18 


3.07 


3.85 




= ^13 




0.28 ± 0.07 


0.34 


0.38 


0.38 


X 


0.44 


0.46 


V\2 






0.1 ± 0.01 


0.12 


0.13 


0.14 


0.2 


0.13 


0.13 
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Having identified the optimal form of the Eshelby matrix, the attention is 
now focused again on the general formulation and the orientation averaging in 
particular, recall Equation (14). Unlike in Section 4.2, however, the joint probability 
density function is represented by real histograms of the fiber tow orientation angles 
already introduced in Figure 1(f); see also Figure 5 for additional examples. With 
such probabilistic characterization in hand, the warp stiffness can be estimated as, 
cf. Equation (15), 

m 

((tf 1 ™*)) = 5>^ (1) (O,0i,O), (22) 

i=i 

where m denotes the number of sampling values and the discrete angles 9i and 
probabilities pi follow directly from the image analysis data. The rest of the analysis 
exactly duplicates the perfect unit cell case. 

The complete data from the analyzed C/C sample involve eleven such histograms 
describing the waviness of the fiber tow in individual plies. The final homogenized 
properties together with the elementary statistical characterization are summarized in 
Table 6. 

Table 6. Homogenized effective stiffnesses determined by the M-T scheme for 
C/C system and histograms measured in (Vopicka 2004). 



Histogram 


rM-T 
^11 


rM-T 
-^12 


rM-T 
-^13 


rM-T 
-^33 


rM-T 
-^44 


rM-T 
-^66 




[GPa] 


[GPa] 


[GPa] 


[GPa] 


[GPa] 


[GPa] 


1 


86.94 


16.71 


15.30 


50.28 


21.63 


42.81 


2 


88.19 


17.16 


15.64 


51.15 


22.23 


44.27 


3 


88.34 


17.14 


15.67 


51.20 


22.24 


44.49 


4 


86.73 


16.64 


15.24 


50.11 


21.49 


42.58 


5 


86.25 


16.47 


15.11 


49.45 


21.27 


42.02 


6 


86.57 


16.58 


15.20 


49.99 


21.63 


42.39 


7 


88.52 


17.28 


15.72 


51.35 


22.35 


44.69 


8 


90.26 


17.90 


16.20 


52.59 


23.23 


46.71 


9 


86.11 


16.42 


15.07 


49.69 


21.21 


41.83 


10 


86.82 


16.67 


15.26 


50.17 


21.54 


42.68 


11 


87.96 


17.07 


15.58 


51.00 


22.15 


43.98 


Average 


87.52 


16.91 


15.45 


50.63 


21.91 


43.49 


Standard deviation 


1.26 


0.44 


0.34 


0.91 


0.60 


1.48 



It becomes clear by associating the results in Table 6 with the corresponding 
micrographs that those segments which have more fibers oriented near the direction 
6 = provide a stiffer response than the others. As an illustration, consider e.g. 
histograms No. 2 and 3 in Figure 5 and corresponding stiffnesses in Table 6. Nevertheless 
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this difference is, due to a rather narrow range of the inclination angles, not too 
important, especially if concerning the approximate character of the M-T method. 

7. Conclusions 

In the present work, an efficient numerical method for the homogenization of plain 
weave composites with both ideal and imperfect tow paths based on the Mori-Tanaka 
method has been proposed. The adopted strategy builds on the matching of results 
of the detailed FEM analysis with the micromechanical model. The most pertinent 
conclusions can be stated as follows: 

i) The simplified method is able to deliver the homogenized parameters with values 
comparable with the detailed finite element simulations. The resulting method 
compares well with independent numerical approaches and available experimental 
data. 

ii) The accuracy of the method depends on the shape of an equivalent ellipsoid, which 
represents both geometrical and mechanical effects, such as inter-tow interactions. 
The parameters of the inclusion follow from a well-defined global optimization 
problem and a FEM-generated training set. Moreover, the optimal shape is robust 
with respect to moderate geometry perturbations. 

iii) The method allows us to directly assess the effects of tow waviness quantified by 
histograms of inclination angles, including the statistical characterization of the 
homogenized stiffnesses. 

The future extension of the method will include the treatment of the intrinsic 
porosity of the C/C composite evident from Figure 1(a). In the framework of multi- 
phase Mori-Tanaka approaches, the porosity can be modeled as an additional phase 
characterized in the simplest case by volume fractions or by three-dimensional computer 
tomography data, see (Piat et al. 2006a, Piat et al. 2006b) for related studies. Such 
work is in progress and will be reported separately. 
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Appendix A. Global optimization algorithm 

The computation scheme is based on the genetic algorithm GRADE (Ibrahimbegovic 
et al. 2004) evaluating the Radial Basis Function Network (RBFN) approximation of 
the objective function, see (Kucerova et al. 2005) for more detailed description. The 
algorithm is briefly described in the flow chart depicted in Figure Al. 

{ Initialize J 



Approximation 
of objective fun. 
by RBFN 



Find optimum 

on RBFN 
using GRADE 



Figure Al. Flow chart of the applied algorithm. 

In particular, instead of directly evaluating the objective function -^(^2,^3) defined 
by Equation (19), the GA evaluates its RBFN approximation. When the optimum of 
the approximation is found, the RBFN is enriched with new neurons according to steps 
described in (Kucerova et al. 2005) and the approximation is refined. At this time 
the real objective function is evaluated at several points. This cycle is repeated until 
the two consecutive solutions differ by less than a certain specified value, set to 10~ 2 . 
Moreover, due to the intrinsic randomness of the algorithm, all reported optimization 
results correspond to the optimum of five independent executions. 
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